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Abstract 


Some important advances took place during the last several years 
in the development of genuinely multidimensional upwind schemes 
for the compressible Euler equations. In particular , a robust, high- 
resolution genuinely multidimensional scheme which can be used 
for any of the flow regimes computations was constructed (sec [1- 
3]). This paper summarizes briefly these developments and outlines 
the fundamental advantages of this approach. 
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1 INTRODUCTION 


The efficiency of the existing steady-state multigrid solvers routinely used for 
the flow problems in engineering practice is still very poor. There clearly exists 
a pressing need for more efficient algorithms. 

In order to obtain a truly efficient steady-state solver, some fundamental issues 
concerning different aspects of an algorithm need to be addressed. The recently 
proposed genuinely multidimensional approach towards the construction of 
the discrete schemes for the compressible flow resolves some of these issues. A 
discussion in this regard is the main subject of this paper. 


1.1 Genuinely multidimensional schemes 


The quest for a genuinely multidimensional upwind scheme began more than 
a decade ago. Initially it was motivated chiefly by the expectation that, once 
such a scheme is developed, it will imitate the physics of the fluid flow more 
accurately than the standard dimension-by-dimension approach. It was, how- 
ever, suggested later in [4], [5] that improving the efficiency of the the steady- 
state solvers may be the most important reason for developing the genuinely 
multidimensional approach. 

One of the main difficulties in the numerical treatment of compressible flow is 
the possible presence of shocks in the solution. It is well known that a scheme 
that is both second order accurate and avoids under- and overshoots (which 
may trigger a nonlinear instability) near discontinuities has to be nonlinear. 
Such a scheme has to incorporate the so-called high-resolution mechanism, 
i.e. a smoothness monitor, that is usually implemented in the form of a flux- 
limiter. Initially, such schemes were developed for the one-dimensional case. 
Then, extending this approach to multidimensions was done on a dimension- 
by-dimension basis. The well known fact, however, is that the Gauss-Scidel re- 
laxation is unstable when applied in conjunction with such schemes [6]. There- 
fore, the standard multigrid solvers have to resort to the defect-correction 
technique or multistage Runge-Kutta relaxation and the efficiency of such 
solvers may be poor. A closer look reveals that the standard high-resolution 
discretizations suffer from the following deficiency: the high-frequency error 
components may be (nearly) invisible to the residuals of the discrete equa- 
tions, i.e. the discrete scheme is (nearly) unstable. In turn this means that 
it may be inherently impossible to construct a good smoother (an important 
ingredient of a multigrid solver) using these discrete schemes. 

A genuinely multidimensional advection scheme was constructed in [4,5]. The 
scheme was named “genuinely multidimensional” since it imitates well the 
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anisotropy of the advection phenomena in two dimensions: artificial dissipation 
is added only along the streamline, while the high-resolution mechanism affects 
significantly the cross-flow direction only. The key feature of this scheme is 
the two-dimensional limiter, i.e. the argument of a limiter-function is the ratio 
of finite differences in two different coordinate directions. The scheme was 
formulated in the control-volume context for Cartesian grids and relied on 
the compact 9-point-box stencil. The fundamental advantage of this approach 
is that the two-dimensional high-resolution mechanism does not damage the 
stability properties of the discretization. 

The so-called “residual distribution” (or “fluctuation-splitting”) schemes for 
scalar advection equation on unstructured triangular grids were presented in 
[7]. It was found later that these schemes have some links to the aforemen- 
tioned genuinely two-dimensional control- volume approach for the advection 
equation. Exploration of these links led to the unification of the two approaches 
and resulted in a scheme that incorporated two-dimensional limiters and was 
formulated for unstructured triangular grids. This scheme (like that presented 
in [4], [5]) can be given a purely algebraic interpretation. However, the task of 
extending these ideas to systems of equations appeared to be a complicated 
one. 

Consider a hyperbolic system of partial differential equations in two dimen- 
sions 


U( + Au x 4- B\iy — 0, 


( 1 ) 


where u is the vector of size N and A , B are N x N matrices. The matrices 
A and B in general do not commute. This means that they cannot be di- 
agonalized simultaneously, i.e. the system cannot be written as N advection 
equations. 

A prolonged effort was to represent locally the physics of compressible flow by 
finite number of simple waves using the local gradients of the solution (in the 
spirit of [8]) with intention to apply a genuinely two-dimensional advection 
scheme to each one of the simple waves. However, the schemes constructed in 
this way for the Euler equations suffered from a severe lack of robustness. 

The breakthrough approach that resulted in a robust genuinely multidimen- 
sional scheme, suitable for the computations of the entire range of the flow 
regimes was presented in [1]. Then it was described in more detail including the 
implications for multigrid and extension to 3D in [2] and [3]. The key idea was 
not to try to apply the multidimensional advection schemes to systems, but 
rather the same strategy that was used to construct the scalar scheme. The al- 
gebraic interpretation of the advection scheme played an important role at this 
point. It was crucial to recognize that a certain linear first order scheme based 
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on standard upwind methodology can be used as a basic building block for the 
hyperbolic systems as well as for the scalar advection. The multidimensional 
high-resolution corrections are then applied in a formal way similarly to the 
scalar case. The resulting scheme for the Euler equations was demonstrated 
to produce a very good quality solution for subsonic, transonic and supersonic 
regimes. The approach was called “genuinely-multidimensional” since it can 
be argued that it leads to a discrete scheme whose artificial dissipation is a 
rotationally-invariant differential operator (in other words the artificial dis- 
sipation operator is independent to a certain extent of the grid direction). It 
was not clear if this particular property is of any direct practical importance. 

The constructed high-resolution scheme for the Euler equations relies on a 
compact stencil. The result of this property is a smaller error in smooth regions 
and better resolution of discontinuities, comparing to the standard dimension- 
by-dimension approach. However, in our view, these are only marginal im- 
provements. The fundamental advantage of this approach is that it leads to 
a scheme that combines high-resolution and good stability properties. It was 
demonstrated in [2,3] that the Collective Gauss-Seidcl relaxation is stable 
when applied directly to the resulting high-resolution discretization of the hy- 
perbolic systems. This results in a very simple, efficient and robust multigrid 
solver for the compressible Euler equations, suitable for the entire range of 
flow regimes. 

Some researchers who were previously pursuing other directions adopted the 
genuinely-multidimensional approach proposed in [1 3] and attributed to it 
a term “Positive Matrix Distribution Schemes”. A modification of the un- 
derlying first order scheme (the system N-scheme) aimed at improving the 
discontinuity resolution was proposed by van der Weide and Dcconinck in [9] . 

It should be mentioned that important steps towards the construction of a 
genuinely multidimensional schemes for the Euler equations were made by 
Colella [10,11], LcVeque [12] and Radvogin [13]. However, the nonlinear high- 
resolution corrections in these schemes rely on one-dimensional limiters, which 
introduces some of the dimension-by-dimension flavor. 

Another very interesting approach was proposed in [14]. A discretization for 
the triangular unstructured grids for the Euler equations was developed. The 
problem was that the scheme was linear, i.c. it did not incorporate any non- 
linear high-resolution mechanism. 


1.2 Multigrid for advection dominated problems 


One of the major reasons for the poor efficiency of the standard flow solvers 
(see [15]) is the fact that for advection dominated problems the coarse grid 


3 



provides only a fraction of the needed correction for certain error components. 
It is well known that the steady Euler equations contain two different factors: 
the advection and the Full-Potential type operators. The latter is either of 
elliptic or hyperbolic type depending on the flow regime (subsonic or super- 
sonic). The difficulty mentioned above can be avoided ([15]) by constructing 
a solver that distinguishes between different factors of the system and treats 
each one appropriately. In the subsonic case, for instance, the advection factor 
can be treated by marching and the elliptic factor by multigrid. The efficiency 
of such an algorithm will be essentially the same as that of the multigrid 
solver for the elliptic part only. Such algorithms are referred to as “essentially 
optimal”. An approach to achieve a separation between the co-factors the 
so-called Distributive Gauss-Seidel relaxation was proposed in [15]. It was 
demonstrated in [16] that using this approach one can obtain the essentially 
optimal multigrid efficiency for a staggered-grid discretization of the incom- 
pressible Navier-Stokes equations. It is interesting to mention that a similar 
observation was made earlier in [17]. 

A related approach was proposed by Ta’asan [18] for the incompressible and 
compressible subsonic Euler equations. The staggered-grid discretization is 
based on the canonical variables formulation (see [19]), that expresses the 
partitioning of the steady Euler equation into elliptic and advection factors. 
The essentially optimal multigrid efficiency was demonstrated using this ap- 
proach for subsonic flow and body-fitted grids. A possible limitation of this 
approach may be that it is not directly gencralizable for the viscous flow. 

Another way towards achieving the optimal multigrid efficiency is based on the 
pressure equation formulation of the Euler equations. We describe it briefly in 
this paper as well. This approach is based on a very old idea and is a general- 
ization [20]. Its main virtue is simplicity. It can also be classified as Weighted 
Gauss-Seidel relaxation [15]. An extensive set of numerical computations us- 
ing this scheme together with more details regarding the implementaion is 
reported in [21]). The limitation of this approach, however, is that it is not 
clear so far if it is generalizablc to viscous compressible flow. 

Following the work of Ta’asan, some researchers attempted to apply the idea of 
partitioning the Euler equations towards the construction of discrete schemes 
(see [22] and [23]). It is well-known that the two-dimensional Euler equations 
in supersonic case can be written as four locally decoupled advection equations 
(sec [24]). This property was used as a basis for applying the advection schemes 
to discretize the system in this case. In subsonic case, however, the distinc- 
tion was made between the advection and the elliptic ( “acoustic subsystem” ) 
partitions. The treatment of transonic flow, however, was problematic since it 
required matching of two different discretizations accross the sonic lines. An- 
other drawback of these approaches was that they are cannot be generalized 
to three dimensions (sec [9]). To conclude, the discrete schemes constructed in 
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this way suffered from a lack of robustness and generality. No optimal multi- 
grid efficiency was demonstrated cither. Finally, some of the researchers who 
previously followed this direction adopted the genuincly-multidimensional ap- 
proach proposed in [13] (see [9]). 


1.3 What this paper is about 


In this paper we first present a brief review of the construction of genuinely 
multidimensional schemes for the scalar advcction and the compressible Euler 
equations. We summarize the basic properties of the discretizations, emphasiz- 
ing those that are unique to this approach and arc of fundamental importance 
for practical purposes. 


The separation of the co-factors can be in general achieved in two ways. One 
way is to cast equations into such a form that the different co-factors can be 
discretized separately. The canonical variables approach by Ta’asan [19] can 
be classified as such. The pressure equation based scheme, presented briefly in 
this paper, belongs to this category as well. The key advantage of this type of 
methods is in their simplicity (this is especially true for the pressure equation 
based scheme). The disadvantage, though, may be in the unsufficient gener- 
ality. Another more general way to achieve the optimal multigrid efficiency 
(see [15]) is to discretize the equations in some primitive form and to apply 
a relaxation of the Distributive Gauss-Seidel type. Such relaxation should be 
designed in such a way that it distinguishes between the different co-factors 
of the system and treats each one of them appropriately. It appears, however, 
that in order to achieve this not any discrete schemes are suitable, but only 
those satisfying a certain condition. As it was mentioned before, it is not clear, 
what are the direct practical implications of the genuine multidimensionality 
(or the rotational invariance of the artificial dissipation) property of the ap- 
proach. However, we present in this paper a heuristc argument suggesting that 
the genuine multidimensionality is closely related to the factorizability prop- 
erty of the discretization. The latter is of fundamental practical importance. 
It is necessary in order to construct a Distributive Gauss-Seidel relaxation 
that will allow to decouple the advcction and Full- Potential co-factors of the 
Euler system and thus to obtain an optimally efficient multigrid solver. This 
approach is very general since it does not require to cast the equations into 
any special form. 
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Fig. 1. Triangle. 

2 MULTIDIMENSIONAL UPWINDING 


In this section we briefly review the construction of the high-resolution gen- 
uinely multidimensional upwind schemes for the scalar advection equation and 
for the Euler system. Consider a general triangular grid covering the domain. 
Assume the discrete unknowns are defined at the grid nodes thus defining a 
linear function and allowing us to evaluate the gradients of the current so- 
lution approximation on each triangle. The approach is to construct discrete 
equations which arc to be solved at each node from the portions of the residu- 
als of the equations computed on the triangles having this node as a common 
vertex. In other words, portions of the residual of the equation evaluated on a 
particular triangle contribute to the construction of the discrete equations at 
the nodes of this triangle. The problem is to find the exact rules for this con- 
struction, so that the resulting discrete equations will have certain desirable 
properties. 


2 A Advection scheme 


We consider triangular clement T, and choose two out of the three faces. Then 
we write the advection equation we wish to solve in the local coordinate system 
aligned with these two faces 


u t + au x + bu y = 0 


( 2 ) 


Without loss of generality we can consider a linear constant coefficient equa- 
tion, since in general non-linear case we can linearize the equation on each 
triangle ([25]). 

We can write the discrete equation at the grid node i in the following form 

^ [ ^ijUj ~ 0, (3) 

]¥=' 
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Definition 1 The discrete scheme (3) is said to be of the positive type ifcij > 
0 for all j. 

Numerical soltuions obtained using a positive scheme satisfy a discrete max- 
imimum principle. This property is useful to ensure that the discrete solution 
will be non-oscillatory near discontinuities. 

We shall outline here the construction of a positive advcction scheme. Residual 
of the equation (2) can be represented as a sum of two portions 

r = r x + r y , (4) 


where 


r 1 = -S T au h x , r y = -S T bu h y . (5) 

Residual of the equation (2) on the triangle T contributes to the construction 
of the discrete equations to be solved at each of the three nodes of the triangle 
according to the following residual distribution formulae 

node 1 <— r x (l — sign(a)) 

node 2 <— r x ( 1 + sign(a)) + r y (l — sign(ft)) (6) 

node 3 <— r y (l + sign (6)) 


It easy to see that this construction results in a positive scheme since for any 
real number z we have the following inequality ± 2(1 ± sign( 2 )) > 0. The 
accuracy of such a scheme, though, is only first order. 

Definition 2 A discrete scheme is called linearity preserving if whenever the 
residual r on the triangle T vanishes, the contributions due to this residual 
lead to a zero update of the solution at each of the three nodes of the triangle. 

A linearity preserving scheme is second order accurate. 

Define the following quantities 

r 1 * = r x + r y 4'(g); r y * = r y + r x A>(q)/q (7) 


where q = —r x /r v . In this paper we assume that is the minmod limiter. Sub- 
stituting r E *,r y * instead of r x , r y into (6) we obtain a high-resolution scheme. 
The important feature here is the two-dimensionality of the limiter, i.e. the fact 
that the argument of the limiter-function is a ratio of numerical derivatives in 
two different coordinate direction ([5], [4]). 
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Using the following limiter identity 


r y ty(q) = — r x ^{q)/q, 


( 8 ) 


it is easy to verify that the constructed nonlinear scheme is indeed both positive 
and linearity preserving. 


2.2 Extension to the Euler system 


Consider a hyperbolic system of partial differential equations 


ip + Au x + Exiy — 0. (9) 

The discrete equation approximating (9) at node % can be written as follows 

n 

Ci,i U, - 53 CijUj = 0, ( 10 ) 

Property of positivity can be formally extended to the system case. 

Definition 3 The discrete scheme (10) is said to be of the positive type if the 
matrices Cij (for all j ) have non-negative eigenvalues. 

It is not clear, however, how to generalize the maximum principle for sys- 
tems. No conclusions can be derived from the positivity property unless the 
additional assumption about the symmetry of the matrices C t J is made. The 
energy stability property of the scheme can be demonstrated in this case. 
However, for the Euler system, this would require the use of the symmetriz- 
ing variables formulation, which is non-conservative. This makes the energy 
stability property too restrictive to be of substantial practical importance for 
the Euler equations. 

It is interesting to note though that the standard high-resolution schemes, if 
carefully implemented, are of the positive type. Therefore, we aimed at con- 
structing a genuinely multidimensional high- resolution scheme ([1 3]), which 
is of the positive type as well. 

Assume that the hyperbolic system of equations (9) is written in the non- 
orthogonal coordinate frame aligned with the two of the faces of triangle T 
(Fig. 1). Residuals of the system on triangle T can be represented as a sum of 
two portions 


R = R x + R y , 


( 11 ) 
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where 

R x = -S T Au h x - 
R y = —S T BUy. 


( 12 ) 


Consider the following residual distribution formula 
node 1 <— R X (I — sign(A)) 

node 2 <— R X (I + sign(A)) + R y (I — sign(Z?)) (13) 

node 3 <— R V (I + sign(£?)) 


Assuming that matrix M has a complete set of real eigenvalues (definition of 
the hyperbolicity of a system) it is easy to sec that matrix ±M(7 ± sign(M)) 
is non- negative definite. This means that the scheme defined by (13) is of the 
positive type by construction (as an upwind scheme is expected to be). 

In order to obtain a positive high-resolution genuinely multidimensional scheme 
for a hyperbolic system we may have first to rewrite system (9) (as it was done 
in [1 3]) in a different set of variables. The auxiliary variables v = (s, u, v, p) 7 
(sec [1 3]) arc a good choice for the Euler system (here s is the entropy, u, v 
- the velocity components and p is the pressure). System (9) rewritten in 
varibales v takes the following form 

v t + Av x + B\y = 0 (14) 


where 


u — Tv 


(15) 


As before, we can compute residual r of system (14) on triangle T and represent 
it as a sum of two portions: 


r = r 1 + r y , 


where 


r* = -S T A\ h x \ 
r y = 


(16) 


(17) 
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Considering the following residual distribution formula 


node 1 <— Tt x (I — sign(.4)) 

node 2 <— T[r x (I + sign(.4)) + r y (I — sign(#))] (18) 

node 3 «— Tr y (I + sign(H)), 


we arrive at a positive first order accurate scheme that is identical to (13). 
Introduce the following quantities 

r f = r i + rfty(qi); rf = rf + r x ^(q t )/q t (19) 

with g, = — rf/rf, and i = 1 , . . . , A r Denote by r 1 * and r y ' vectors whose 
components are rf and rf (i = 1, . . . , N) respectively. The high- resolution 
genuinely two-dimensional scheme can be obtained by substituting r x * and r v * 
instead of the r 1 ' and r y into (18). 

Using the limiter identity 


r l V(Qi) = -U I 'J KQi)/Qi> for * = 1, N (20) 


it is possible to show that the genuinely two-dimensional high-resolution scheme 
is both positive and linearity preserving. We emphasize here again, that in 
order to achieve this property for the Euler equations, it was necessary to 
use another (non-conservative) form of the equations when introducing the 
high-resolution mechanism ([1 3]). The auxiliary variables formulation was 
suitable for this purpose. The only justification for the desirability of the pos- 
itivity property for systems of equations is that the standard high-resolution 
schemes, if properly implemented, have this property. 

Remark 4 Recall that when constructing a discrete scheme we have chosen 
arbitrarily two out of three faces of the trinagle T ( see Fig.l) for the purpose of 
numerical derivatives evaluation. Note, that any choice will result in a scheme 
with good properties. In the case of the scalar advection, though, if we use 
two faces which are both of the same kind with respect to the flow direction 
(inflow or outflow), the resulting scheme will rely on the narrowest possible 
stencil. Therefore, it will have a smaller cross-stream error coefficient and will 
provide a somewhat better disontinuity resolution. It is possible in the case of 
the Euler system to optimize the resolution of a specific sharp layer (shock, 
contact discontinuity or shear layer) by the appropriate choice of the faces 
(two “ closest ” to the layer direction. In general, though, it seems reasonable 
to choose two faces which are either inflow or outflow. This will provide better 
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resolution of the contact discontinuities and shear layer, while shocks are well 
resolved anyway. 

In addition to the smaller error in smooth regions and sharper resolution of 
discontinuities (compared to the standard dimension- by-dimension methods), 
the constructed scheme also offers a possibility to optimize the stencil in order 
to resolve a particular discontinuity layer (by choosing two faces of a triangular 
element that arc the closest to the direction of this layer). The important 
advantage of this approach, however, is that the genuinely multidimensional 
approach presented here leads to a scheme which has both high-resolution and 
good stability properties. 


3 SEPARATION OF CO-FACTORS 


Incompressible steady Euler equations in two dimensions can be written in 
the following matrix form (assuming that the fluid density p = 1) 

Lu = 0, (21) 


where 


u = (u,v,p) T , 


‘ Q 0 ^ 
0 Q d y 
yd x d y 0 j 


( 22 ) 


(23) 


and Q = U • V is the advection operator. In order to find out what is the type 
the system of equations we can look at the determinant of matrix L 

det(L) = -Q 2 A (24) 


The co-factors of the determinant are of the advection and another of the 
elliptic types. 

Compressible steady Euler equations in two dimensions can be written in the 
following matrix form 


Lu = 0, 


(25) 
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where 


u = ( S,U,V,p ) T , 

1 Q 0 0 0 ^ 

0 pQ 0 d x 

0 0 pQ d y 

\ 0 pdx pd y £ Q ) 


(26) 


(27) 


The determinant of matrix L in this case 


det(L) = p 2 (Q) 2 [i<3 2 - A] 
c z 


(28) 


The determinant has two distinct co-factors: one of advection and another of 
the Full- Potential type. 

In the rest of this paper we shall consider the subsonic regime only. The Full- 
Potential factor in this case is of the elliptic type. 

It was suggested in [15] that the different co- factors should also be treated 
differently (each one in the appropriate way). One way to do this is to cast the 
equations into such a form that the co-factors can be discretized separately. 
Canonical variables formulation [18] as well as the pressure equation based 
schemes, that are described in the next section, belong to this category. A more 
general approach (as proposed by Brandt in [15]) is to discretize the equations 
in some primitive form, but to design relaxation of Distributive Gauss-Seidel 
type (DGS) such that it will separate the treatment of co-factors. However, the 
discrete scheme suitable for this purpose should be factorizable, i.c. satisfy the 
discrete analog of the property (28). We show in this paper that the genuinely 
multidimensional upwind approach leads to a factorizable scheme. 


4 PRESSURE EQUATION APPROACH 


In this section we briefly describe the pressure equation based approach first 
for incompressible and then for compressible cases. 
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4-1 Incompressible case 


Considering a triangular (unstructured) grid and assuming that the unknowns 
are located at grid vertices, we can define pieccwize-lincar functions approxi- 
mating each unknown. Residual of the Euler equations can then be evaluated 
on each triangular element of the grid: 


r = L h n h , 


where = (u h ,v h ,p h ) T , 


L h 


( Q h 0 d h \ 
0 Q h 

d v 0 j 


(29) 


(30) 


and 


Q h = u h d £ + v h dy 


(31) 


is the discrete advcction operator. 

We would like then to construct the discrete equations to be solved at a certain 
grid node by assembling in a certain way the residuals of the Euler equations 
computed on the triangular element having this node as a common vertex. 
These equations can be written in the following form 

P h r = P h L h u h - 0 (32) 


In order to obtain the equations with desirable properties, i.e. to achieve the 
decoupling of the co-factors, we can define the following assembly matrix P h 


ph 


0/0 

m h ) 


(33) 


where 


mHf h ) = (u h f ) - aft®*/*) 04) 
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is the operator adjoint to Q h . Then 



{ Q h 0 <9^ 


( h\ 
U 

P h r = 

0 Q h d h y 


yk 


K 0 0 A h ) 




+ s.p.t. 


(35) 


where A h is the discrete Laplacian and “s.p.t.” stands for subprincipal terms 
which appear due to the non-constancy of the advection coefficients and are 
not important for the purpose of constructing a relaxation scheme [15]. The 
maxtrix of the finite difference operators in (35) is upper-triangular. The pres- 
sure is subject to a Poisson equation wich is decoupled (upto subprincipal 
terms) from the rest of the system. The standard Gauss-Seidel relaxation 
can, therefore, be applied to it as a smoother. The momentum equations can 
be looked at as advection equations with (known) forcing terms (pressure 
derivatives). The ideas of multidimensional upwinding can be applied to dis- 
cretize them. The strategy applied to relax the system can be then to relax 
the pressure first and then to update the velocities components by relaxing 
the momentum equations in the downstream direction. An extensive set of 
computational results using this approach is presented in [21]. 


4-2 Compressible case 


Similarly to the incompressible case, residuals of the equations can be evalu- 
ated on each triangular element: 


r = L h u\ 


where u h = (s h , u h , v h ,p h ) T and 

( Q h 0 0 0 

h = 0 p h Q h 0 d h x 

0 0 p h Q h d* 

V 0 P h d* P h dy cbQ h ) 


(36) 


(37) 


and & is the (discrete) speed of sound. Again, we construct the discrete equa- 
tions to be solved at each node assembling in a certain way the residuals of 
the equations on the elements surrounding this node 

ph r = ph L h n h = q (3g) 
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Choosing the assembly matrix to be 


ph 


f I 0 0 0 ^ 

0/0 0 
0 0 1 0 
VO %d£(Q*) A J 


(39) 


we obtain 


P h r 


1 Q h 0 0 

0 p h Q h 0 
0 0 p h Q h 

^00 0 



+ s.p.t. 


(40) 


The principal part of last equation is the discrete Prandtl-Glauert (or Full 
Potential) operator acting on the pressure. This operator is of the elliptic type 
in the subsonic regime. The solution process of the resulting discrete equa- 
tions is very similar to that for the incompressible case. Some numerical tests 
illustrating the efficiency of the multigrid solver based on this discretization 
arc presented in the §6.2. 


5 UPWINDING AND CO-FACTORS SEPARATION 


Now we return to the genuinely multidimensional upwind approach. We would 
like to construct a linear first order upwind “positive” scheme such that it is 
factorizablc and is also upgradable to second order using the genuinely two- 
dimensional high-resolution mechanism. 

First, we shall take a closer look at the one-dimensional case. 


5.1 One- dimensional case 


Consider a first order upwind scheme for the one-dimensional Euler equations. 
Without loss of generality we consider the primitive variable formulation 

L h u h = 0 (41) 
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where 


u h = (s h , u h , p h ) T (42) 

and 

%\u\d* x + Q 2h 0 0 ^ 

0 p{-\ cd h xx + Q 2h ) d 2 x h -^ x , (43) 

0 p(af-|^) -A dxx + ^Q2 hj 

where h is a meshsize, d xx is a central approximation of the second derivative, 
d 2h is a central approximation of the first derivative and Q 2h = ud 2h is the 
advection operator. 

5.1.1 Factorization 
The determinant of L h : 

iet(L h ) = \u\dt, + Q 2 ‘)[(l - M 2 )<&] (44) 

The first factor is the upwind scheme approximating an advection operator 
corresponding to the entropy equations. The Full-Potential factor is approxi- 
mated by a “short” central difference. The issue of factorization appears to be 
trivial in this case, since the momentum and the pressure equations correspond 
solely to the elliptic factor. 

5.1.2 Distributive Gauss-Seidel relaxation 
Introducing new variables 

(*\ u h , p h ) T = u h = M h w h = M h (s h , w h , (f) h ) T (45) 

the Euler system will take the following form 

L h M h w h = 0 (46) 
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Assuming that that 


M h = 


1 1 0 0 
{0-p(dl l '-l*dZ,)p(ic% x -Q™)) 


(47) 


we obtain 


L h M h = 


( -\\ u \ &xx + Q 2h 

o 

o 


o o ^ 

p{ 1 - M 2 )% x 0 

0 p(l - M*)d h xx ) 


(48) 


The philosophy of the Distributive Gauss-Seidcl relaxation is as follows: We 
would like to store at the gridpoints the primitive variables u and not the 
auxiliary variables w. For this purpose the updates in w and cf> corresponding 
to the relaxing the elliptic factor at point i should be translated into the 
updates in u,p at points ( i — 1), (i), ( i + 1) according to M h . 


5.1.3 DGS relaxation and the Riemann solver 

It is well known that the one-dimensional Euler system can be diagonalized, 
i.c. rewritten as a set of (locally) decoupled advection equations for the char- 
acteristic variables (s,a + ,a~) r , where 

a x — pcu x + p x , and a~ = pcu x - p x . (49) 


Some algebra reveals that relaxing the Full-Potential (elliptic) factor corre- 
sponds to: 

at point (i — 1) - update o: + , keep a - the same; 
at point (i + 1) - update a - , keep a + the same; 
at point i - update both a + and a~. 

We can conclude that there arc some links between the characteristic variables 
formulation (approximate Riemann solver) and the design of the DGS. 
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5.2 Two-dimensional case 


We shall look now for a factorizable upwind scheme in two dimensions. 


5.2.1 Dimension-by- dimension approach 

Considering here the isentropic case (no loss of generality for the purpose of 
the discussion presented here), we can write such a scheme in the following 
symbolic form 


L h u h = 0 


(50) 


The modified equations (or the First Differential Approximation FDA) cor- 
responding to the scheme dimension-by-dimension 


FDA{L h ) = 

' p [- 1 (cdxx + \v\dyy) + Q] 0 d x 

0 p[-%(\u\d xx + cdyy) + Q\ dy 


| \ud xx 


— — Ijf) 

2 c u yy 


(51) 


V P(d x -^ c d xx ) p { dy -bvQ yy) _|I A + Q ) 

It is easy to sec that the matrix (51) cannot be factorized. 


5.2.2 Genuinely multidimensional approach 

The approach towards the construction of discrete schemes for the Euler equa- 
tions ([1,2]) was called “genuinely multidimensional” since it leads to schemes 
that retain (to a certain extent) the rotational invariance property of the Eu- 
ler equations. Namely, it can be argued that the artificial dissipation terms 
present in these schemes approximate a rotationally invariant differential oper- 
ator. In its turn this may mean that the waves oblique to the grid are “properly 
upwinded” or, in other words, the same approximate Ricmann solver can be 
“recovered” in an arbitrary direction. 

It is not clear whether or not this property, though intuitively appealing, is 
of any direct practical importance for the steady-state computations. There- 
fore, we do not discuss it in detail. However, we have observed previously that 
there are some links between the approximate Riemann solver and the design 
of DGS in one dimensional case. Therefore, the following conjecture seems 
reasonable. 
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Conjecture A genuinely multidimensional scheme is factonzable. 


It was pointed out in [2] that some of the multidimensional second order correc- 
tions can be added without limiters resulting in a linear ‘positive’ scheme with 
essential multidimensional character. Namely, for the w-rnomentum equation 
in subsonic case, those are the cross- derivative correction terms that compen- 
sate for the loss of accuracy due to the artificial dissipation in x direction. For 
the u-momentum equations those will be the terms that compensate for the 
loss of accuracy due to the artificial dissipation in y-direction. 

Writing such a scheme in the symbolic form 

L[ 2D] u h = 0, (52) 


and considering the corresponding FDA 


FDA{L\ 2D] ) = 

f p[-\{cdxx + \v\dyy) + Q) p[-%{c- \v\)d X y\ d x - \ \Qx ^ 

p\-\{C - \u\)d xy \ p[-%(\u\d xx +cdyy) +Q\ dy-^ c Qy 

V p{d x -^ c d xx ) p(dy -hv cdyy) _|i A + q) 


(53) 


we can easily verified that FDA matrix is factorizable. The added multidimen- 
sional correction played a crucial role in achieving this property. 


However, the FDA does not uniquely define a discrete scheme. Moreover, not 
all the discrete schemes corresponding to a certain FDA are factorizable. A 
factorizable scheme corresponding the above mentioned FDA was constructed 
on a quad-type grid. The details concerning the scheme as well as the Dis- 
tributive Gauss-Scidel relaxation will be given elsewhere. 


6 NUMERICAL EXPERIMENTS 

6.1 Multidimensional upwinding 


The purpose of the numerical experiments reported in this section is to demon- 
strate the robustness of the genuinely multidimensional upwind scheme and 
the quality of the numerical solutions obtained by its means. The multigrid 
algorithm employs the lexicographic Gauss-Scidel relaxation. 
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a) 



b) 



Fig. 3. Supersonic flow in a channel over a circular bump, grid 161 x 33: a) solution 
obtained by 2 FMF — H’(2. 1) algorithm; b) the same, except that 3 more W{ 2, 1) 
cycles were performed on the finest grid. 


It was demonstrated in [4], [5] that the IF MG — W (2, 1) algorithm employing 
the genuinely two-dimensional advection scheme (and Gauss-Seidel relaxation 
with direction-free ordering - like Red-Black) is capable of producing second 
order accurate (both in smooth regions and in terms of discontinuity location) 
solution to such a problem. We expect this to be true for the Euler equations 
as well (though more studies should be performed). Therefore, we present 
solutions obtained using this algorithm for a few tcstcases. 


6.1.1 Supersonic flow in a channel with a bump. 

The test case considered here is a supersonic (Mach=2.9) flow in channel with 
a circular bump. The bump is located at the lower wall of the channel at 
1 < x < 2 and its surface is a circular arch of 7r/3 and radius 1. Note that 
the actual shape of the domain is a rectangle. The influence of the bump on 
the flow is imposed through the boundary conditions: the velocity component 
normal to the surface of the bump at a certain location is being reflected. 

The experiment uses the finest grid of the size 161 x 33 points. The solution 
obtained by 2 FMG — W( 2,1) algorithm is presented on Fig.3(a). Fig. 3(b) 
presents the numerical solution obtained using the same algorithm but per- 
forming 3 more cycles (total 5) on the finest level. 
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Fig. 4. Transonic flow over a wall with a circular bump: free stream Mach=.9, grid 
200 x 200 pts. 

6.1.2 Transonic flow over a circular bump 

The testcasc considered here is concerned with a transonic flow over a flat 
wall with a bump. The surface of the bump is again a circular arch of 7r/3 and 
radius 1 and its location is between 3.5 < x < 4.5. Again, in order to keep the 
experiments simple at this stage of work, the bump is treated the same way 
as in the previous experiments. The free stream Mach=.9 in this case. The 
shock of the “fish-tail” shape is generated in this case (Fig. 4). 


6.1.3 Subcritical flow past an airfoil. 

5 Here we present an experiment concerned with the subcritical flow past an 
airfoil. The testcasc considered is Mach = .63 flow past AMCA0012 airfoil at 
the angle of attack of 2°. The grid conatins about 9800 nodes. Pressure and 
density contours arc presented by Figs. 5 (a) and (b) respectively. 


6.2 Pressure equation based schemes 


Here we illustrate the efficiency of the multigrid algorithm that distinguishes 
between the different co-factors of the system. The testcasc is a flow in a 
channel with a bump. The mutligrid cycle employed is V' (2, 1), the relaxation 
is lexicographic Gauss-Seidel in the flow direction. The solution contour plots 
for the incompressible case are presented on Fig.6. The sample convergence 
rates for the incompressible and compressible subsonic cases can be found in 
Fig.7. It can be obserevd that the residuals reduction per multigrid cycle is 
almost an order of magnitude. Also, the extensive set of computational results 
presented in [21] indicates that this behavior is essentially grid-independent. 
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b) 



Fig. 5. Subcritical flow (Mach = .63) over NACA0012 airfoil, angle of attack - 2°, 
grid of ~ 9800 nodes: a) pressure contours; b) density contours. 

7 CONCLUSIONS 


An approach towards constructing a genuinely multidimensional upwind scheme 
was intoduccd in [1 3] . The fundamental advantage of this approach for prac- 
tical purposes is that the multidimensional high-resolution mechanism (unlike 
the standard one) does not damage the stability properties of the scheme. The 
conclusion made in this paper is that the genuinely multidimensional approach 
leads to a scheme that is also factorizable. The practical importance of this 
property is that an optimally efficient multigrid solver can be obtained through 
the construction of an appropriate Distributive Gauss-Seidel relaxation, that 
distinguishes between the different co- factors of the equations. Also, since the 
factorizability property can be easily verified, we suggest it is used as the 
definition of genuine multidimensionality of a scheme. 
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b) 



Fig. 6. Incompressible flow in a channel with a bump, grid of 97 x 33 vertices: a) 
v-vclocity contours; b) pressure contours. 



Fig. 7. Sample convergence histories for incompressible and compressible (Mach=0.2 
and Mach=0.5) cases. 

Wc also presented briefly the pressure equation based schemes and demon- 
strated on their example what efficiency of the multigrid solver can be ex- 
pected if the distinction between the different co-factors of the equations is 
made by the algorithm. 

Due to its generality, the genuinely multidimensional approach for discretiza- 
tion of the Euler equations may play a crucial role in constructing a general 
optimally efficient multigrid flow solver suitable for engineering computations. 
This is because 

it does not rely on casting the equations into any special form; 
extends to the compressible Navier-Stokes equations. 
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